package scu.maqiang.fes;

import static scu.maqiang.numeric.MVO.average;
import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES3T4PE extends FES {

	public FES3T4PE(Mesh3T4 mesh) {
		super(mesh, 4);
		eleStiff = FES3T4PE::elementPEStiff;
		eleMass = FES3T4PE::elementPEMass;
		eleSource = FES3T4PE::elementPESource;
		eleFlux = FES3T4PE::elemenSurfaceForce;
	}

	/**
	 * 计算单元压电刚阵, 注意排列顺序为11, 22, 33, 23, 31, 12, 1, 2, 3
	 * 
	 * @param xyz       单元坐标数组
	 * @param mat       5x5 本构矩阵，下标排列顺序为11, 22, 33, 23, 31, 12, 1, 2, 3
	 * @param tp        问题类型
	 * @param K 得到的单元刚度矩阵
	 */
	public static void elementPEStiff(int[] ele, double[][] xyz, double[][] mat, BVPType tp, double[][] K) {
		double[][] grad = new double[4][3];
		double vol = CG.gradient3DT4(xyz, grad);
		double[][] B = null;
		MVO.fill(K, 0.0);
		switch (tp) {
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			B = new double[9][16];
			for (int i = 0; i < 4; i++) {
				B[0][4 * i] = B[5][4 * i + 1] = B[4][4 * i + 2] = B[6][4 * i + 3] = grad[i][0];
				B[1][4 * i + 1] = B[5][4 * i] = B[3][4 * i + 2] = B[7][4 * i + 3] = grad[i][1];
				B[2][4 * i + 2] = B[3][4 * i + 1] = B[4][4 * i] = B[8][4 * i + 3] = grad[i][2];
			}
			mul(matmul(transpose(B), matmul(mat, B)), vol, K);
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assemblePEMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES3T4PE::elementPEMass, tp, A);
	}

	public void assemblePEMass(double[] constCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(constCoef, FES3T4PE::elementPEMass, tp, A);
	}

	public static void elementPEMass(int[] ele, double[][] xyz, double[][] r, BVPType tp, double[][] K) {
		double[] rho = r[0];
		double sumRho = MVO.sum(rho);
		double volume = CG.volume3DT4(xyz);
		MVO.fill(K, 0.0);
		if (tp == BVPType.CONSIST_MASS) {
			K[0][0] = K[1][1] = K[2][2] = (2 * rho[0] + sumRho) / 60 * volume; // 11
			K[4][4] = K[5][5] = K[6][6] = (2 * rho[1] + sumRho) / 60 * volume; // 22
			K[8][8] = K[9][9] = K[10][10] = (2 * rho[2] + sumRho) / 60 * volume; // 33
			K[12][12] = K[13][13] = K[14][14] = (2 * rho[3] + sumRho) / 60 * volume; // 44

			K[0][4] = K[1][5] = K[2][6] = K[4][0] = K[5][1] = K[6][2] = (rho[0] + rho[1] + sumRho) / 120 * volume; // 12
			K[0][8] = K[1][9] = K[2][10] = K[8][0] = K[9][1] = K[10][2] = (rho[0] + rho[2] + sumRho) / 120 * volume; // 13
			K[0][12] = K[1][13] = K[2][14] = K[12][0] = K[13][1] = K[14][2] = (rho[0] + rho[3] + sumRho) / 120 * volume; // 14
			K[4][8] = K[5][9] = K[6][10] = K[8][4] = K[9][5] = K[10][6] = (rho[1] + rho[2] + sumRho) / 120 * volume; // 23
			K[4][12] = K[5][13] = K[6][14] = K[12][4] = K[13][5] = K[14][6] = (rho[1] + rho[3] + sumRho) / 120 * volume; // 24
			K[8][12] = K[9][13] = K[10][14] = K[12][8] = K[13][9] = K[14][10] = (rho[2] + rho[3] + sumRho) / 120
					* volume; // 34

		} else {
			K[0][0] = K[1][1] = K[2][2] = (rho[0] + sumRho) / 20 * volume;
			K[4][4] = K[5][5] = K[6][6] = (rho[1] + sumRho) / 20 * volume;
			K[8][8] = K[9][9] = K[10][10] = (rho[2] + sumRho) / 20 * volume;
			K[12][12] = K[13][13] = K[14][14] = (rho[3] + sumRho) / 20 * volume;
		}
	}

	public static void elementPESource(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double volume = CG.volume3DT4(xyz);
		double[] fu = f[0];
		double[] fv = f[1];
		double[] fw = f[2];
		double[] rf = f[3];
		double sumFu = MVO.sum(f[0]);
		double sumFv = MVO.sum(f[1]);
		double sumFw = MVO.sum(f[2]);
		double sumRf = MVO.sum(f[3]);
		Arrays.fill(eleF, 0.0);
		if (tp == BVPType.COMMON) {
			for (int i = 0; i < 4; i++) {
				eleF[4 * i] = (fu[i] + sumFu) * volume / 20.0;
				eleF[4 * i + 1] = (fv[i] + sumFv) * volume / 20.0;
				eleF[4 * i + 2] = (fw[i] + sumFw) * volume / 20.0;
				eleF[4 * i + 3] = (rf[i] + sumRf) * volume / 20.0;
			}
		}
	}

	public static void elemenSurfaceForce(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] fu = f[0];
		double[] fv = f[1];
		double[] fw = f[2];
		double volume = CG.area3DT3(xyz);
		double sumFu = MVO.sum(fu);
		double sumFv = MVO.sum(fv);
		double sumFw = MVO.sum(fw);
		for (int i = 0; i < 3; i++) {
			eleF[4 * i] = (sumFu + fu[i]) * volume / 12;
			eleF[4 * i + 1] = (sumFv + fv[i]) * volume / 12;
			eleF[4 * i + 2] = (sumFw + fw[i]) * volume / 12;
		}
	}

	public void assembleSurfaceElectricDisp(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES3T4PE::elemenElectricDisp, tp, RHS, bdLabel);
	}

	public void assembleSurfaceElectricDisp(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES3T4PE::elemenElectricDisp, tp, RHS, bdLabel);
	}

	public void assembleSurfaceElectricDisp(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES3T4PE::elemenElectricDisp, tp, RHS, bdLabel);
	}

	public static void elemenElectricDisp(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] eleF) {
		double[] Rf = g[0];
		double volume = CG.area3DT3(coord);
		double sumRf = MVO.sum(Rf);
		for (int i = 0; i < 3; i++) {
			eleF[4 * i + 3] = (sumRf + Rf[i]) * volume / 12;
		}
	}

}
